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Abstract: Sterile Neutrinos with a mass in the keV range form a good candidate for 
dark matter. They are naturally produced from neutrino oscillations via their mixing with 
the active neutrinos. However the production via non-resonant neutrino oscillations has 
recently been ruled out. The alternative production via Higgs decay is negligibly small 
compared to neutrino oscillations. We show that in the neutrino-phillic two Higgs doublet 
model, the contribution from Higgs decay can dominate over the contribution from neutrino 
oscillations and evade all constraints. We also study the free-streaming horizon and find 
that a sterile neutrino mass in the range of 4 to 53 keV leads to warm dark matter. 
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1 Introduction 

The Standard Model (SM) of particle physics is very successful, but it fails to explain 
neutrino mass and dark matter (DM). Dark matter accounts for about one quarter of 
the energy density of the Universe, five times more than ordinary matter, but its origin is 
unknown. A good candidate for the dark matter are sterile neutrinos with a keV-scale mass 
and tiny mixing with the active neutrinos, which is a simple extension of the SM [1, 2]. 
In contrast to standard cold dark matter, they are generally warmer with a larger free- 
streaming horizon. Thus they are candidates for warm dark matter and suppress structure 
at small scales addressing the missing satellite problem [3-5] and possibly explaining the 
velocities of pulsars [6, 7]. 

There are many ways of producing sterile neutrinos: (i) they can be produced through 
neutrinos oscillations in the early Universe via a small mixing with the active neutrinos [8- 
10]. This mechanism is already excluded by observation [11], but the bounds can be 
avoided, if there is a large enough primordial lepton asymmetry and sterile neutrinos are 
produced via resonant oscillations [12]. (ii) Another well-studied alternative is non-thermal 
production via decay of a scalar field in thermal equilibrium [6, 13-17], or a scalar produced 
via the freeze-in mechanism [18], which subsequently decays to sterile neutrinos [17, 19-21]. 
Recently several alternative production mechanisms from decay have been suggested such 
as the production from the decay of pions [22] , Dirac fermions [23] , light vector bosons [24] , 


- 1 - 




or a condensate formed during inflation [25, 26]. (iii) Finally the keV sterile neutrinos 
could have been In thermal equilibrium and their abundance diluted by production of 
entropy [27-31]. 

In any model with mixing between active neutrinos and the keV sterile neutrinos 
a fraction of the sterile neutrino abundance will be generated via neutrino oscillations. 
The mixing is generally induced after an electroweak doublet scalar obtains a vacuum 
expectation value (vev). The Yukawa interaction however induces a second contribution: 
The electroweak doublet scalar can decay into a SM lepton and a sterile neutrino. In 
a model with one Higgs doublet, this contribution is always negligible compared to the 
contribution from neutrino oscillations [32], because the vev, v = 174 GeV, and thus the 
mixing is sizeable. However this does not hold anymore in models with multiple Higgs 
doublets. The vev of one of the Higgs doublets might be tiny, smaller than C)(MeV), and 
consequently the production via Higgs decay might dominate. Sterile neutrino dark matter 
with a keV-scale mass has been previously considered in a two Higgs doublet model in 
Ref. [33] and the production of the required number density via the decay of an electroweak 
doublet has been studied in Ref. [34] in the scotogenic model of neutrino mass [35]. 

We consider a two Higgs doublet model, where one of the electroweak doublet scalars 
exclusively couples to the sterile neutrino and study the production of keV sterile neutrino 
DM via the decay of this electroweak doublet in detail. The main result is the momentum 
distribution function for the sterile neutrino DM and the free-streaming horizon, which 
we use to determine the relevant parameter space where the keV sterile neutrino consti¬ 
tutes warm dark matter. This mechanism can be easily embedded in a seesaw [33, 36] or 
radiative [34, 35] neutrino mass model. 

The paper is organised as follows: In section 2, we introduce the two-Higgs doublet 
model and discuss the mass spectrum. The produced sterile neutrino DM abundance is 
discussed in section 3. Section 4 is dedicated to the free-streaming horizon of the sterile 
neutrinos and we briefly comment on the effective number of relativistic degrees of freedom 
in section 5. In section 6, we discuss the possibility that the scalar doublet obtains a tiny 
vev, thus the sterile neutrinos can decay and explain the observed an X-ray line at 3.55 
keV [37, 38]. We conclude in section 7. Technical details are collected in the appendices. 

2 Neutrino-phillic Two Higgs Doublet Model 

We consider a two Higgs doublet model with a second scalar doublet Hi, which exclusively 
couples to the sterile neutrino N and the left-handed lepton doublet L. This is guaranteed 
by a Z 2 symmetry under which the new fields, N and Hi^, are odd, but all SM particles are 
even. A coupling to other fermions is strongly constrained by flavour-violating processes. 
The most general Yukawa interactions in the lepton sector are given by 

- £ = tjeLHE^ + ULNLHyN + + h.c. . (2.1) 
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and the most general scalar potential is defined in the usual way 


y = (2.2) 

+ h.c.] , 

where we used a field redefinition of H^, to absorb the complex phase of A 5 . The Yukawa 
couplings yLN,a are generally complex. After the SM Higgs doublet obtains a vev, = 
mf/Ai, electroweak symmetry is broken and we decompose the fields in terms of their 
components 

" “ I’' + 7! ('> + «“)) “ (a ^ ^ 

The scalar masses at leading order are given by 

ml = 2m\ = 2Aiu^ = m\ + (A 3 + A 4 + A 5 ) (2.4) 

m]^o =ml + (A 3 + A 4 - As) m|-± =ml + , 

where h describes the observed Higgs boson at mh = 125 GeV [39, 40]. As long as Hi, 
does not obtain a vev, there is no mixing between the different states. We will comment 
on this possibility in section 6 . The active neutrinos obtain mass in the usual way. Both 
Dirac and Majorana mass terms are possible. For example, one can introduce the Weinberg 
operator[41], LLHH, to generate the neutrino mass. We will not discuss it further, because 
it does not affect the production of the keV sterile neutrino N. 


3 Dark Matter Production 


Freeze-in production of sterile neutrino DM with a second Higgs doublet has been studied 
in the scotogenic model [34] . Here we focus on keV sterile neutrinos and do not specify the 
mechanism of neutrino mass generation explicitly. 

The keV sterile neutrino can be produced via the decay of the scalar fields (A:, K^, H^) 
while they are in thermal equilibrium. We will present a crude simple calculation using 
the Maxwell-Boltzmann approximation and neglect Pauli-blocking. The result will result 
in a good order of magnitude estimate. In our discussion, however, we will use the more 
accurate result using the distribution function, which can be found in App. C. Assuming 
that inverse decays can be neglected, the yield Y{T) = n{T)/s{T) of the sterile neutrino 
can be obtained from the Boltzmann equation. 


dYj^ ^ 7Afi (T) 

® dT H{T) ’ 


(3.1) 


where s is the entropy density of the Universe, H(T) is the Hubble parameter at a given 
temperature and 7 Ar^(T) is the thermally averaged sterile neutrino production rate. 


7Wi(T) = ^ N,l), 

X 


(3.2) 
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where X = k, K^, and I is a SM lepton. Following the freeze-in calculation in Ref. [18], 
we can integrate the equation and obtain the final yield of sterile neutrinos after freeze-in 


woo _ 45 


X 


gxTjX ^ Nil)Mo 
m\gl{Td) 


x^Ki{x)dx 


where we defined the typical mass scale 


(3.3) 


Mo 


3 



(3.4) 


This allows us to rewrite the Friedmann equation during the radiation dominated epoch as 


H = — 
Mo 


(3.5) 


using the Planck mass Mpi and the usual definition of the effective entropy degrees of 
freedom g^ and the relativistic degrees of freedom g'^ 


s 


2 ^ 

45 






(3.6) 


We can obtain an approximate analytic solution to this equation by extending the integra¬ 
tion over all positive values of x, i.e. Xmin —^ 0 and Xmax —^ oo 


iq^ 9/2 ^ my^(Td,x)Vg*(W 

(3.7) 

0.328 V 

X '^x9*i'^d,x)V9*{Td,x) 

(3.8) 


In order to obtain simple analytic results, we made several (crude) approximations: (i) 
We used the Maxwell-Boltzmann approximation and thus also neglected Pauli-blocking of 
the neutrinos, (ii) We extended the integration boundaries in (3.3) to obtain the leading 
order result in Eq. (3.7). This is justified by noting that freeze-in is typically dominated 
by processes around T ~ mx[18]. (iii) We assumed that the effective number of relativistic 
degrees of freedom for entropy, g^, and energy, g^, do not change during the production of 
dark matter, which is reasonably well satisfied for scalar masses mx ^ 100 GeV. (iv) We 
neglected 2 0 2 scattering processes like X^ + ^ N + and X^ + ^ N + . 

These processes are subdominant compared to two body decays of X due to phase space 
suppression, (v) We assumed that the particles X are in thermal equilibrium until freeze-in 
occurs at x/i ~ 2 — 5 [18] and thus Tfi > 20 GeV for mx ^ 100 GeV. For sufficiently large 
Higgs portal couplings^ Aa^s, the scattering of the scalars X with b quarks, X + b ^ X + b, 
will keep the scalars X in thermal equilibrium similar to the SM Higgs down to T ~ 5 
GeV. (vi) Finally we use the usual vacuum decay rate and do not take into account finite 
temperature effects, which has been studied e.g. in Ref. [42]. We expect these corrections 

^Note that the Higgs portal couplings enter the freeze-in calculation only indirectly via the mass mx of 
the particle X. 
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(a) mjv vs. -x/X^ct \yLN,a\^- The contour labels 
are the scalar mass irikk in GeV. 



rritf [keV] 

(b) rrikk vs. 171 ^. The contour labels are the 
effective coupling \/Xq \yLN,a\‘^- 


Figure 1: Contour plots with fixed DM relic abundance = 0.1199 [43]. 


to be small, because the Yukawa coupling is small and the sterile neutrino a gauge singlet. 
In the following we will, however, use the more accurate result in Eq. (C.l). 


Voo 


~ 0.207^ 

X 


gxTjX ^ Nil)Mpi 

{Td,x ) V g* {Td,x ) 


(3.9) 


It has been obtained by solving the Boltzmann equation for the distribution function 
without using the Maxwell-Boltzmann approximation and including Pauli-blocking. The 
detailed calculation is outlined in App. B and C. 

The number of degree of freedom of the scalar fields are given hy gk = 1, gxo = 1, gK+ = 
2 and the decay widths T^ are given in App. A. We finally obtain the relic abundance of 
the sterile neutrino dark matter using /pcr with the critical energy density 

Per = 3H‘^Mpi/8TT and find 


^ 6.88 X 


lo^^mAT 

X 


gxT{X ^ Nil) 

m^gl{Td,x)\/g* {Td,x) 


(3.10) 


Taking the limit of equal scalar masses, ~ m^o — 'mK+ = T^ x = T^, the relic 

abundance 




6.88 X 10^® Smx , 

gi{Td)^/g^*{Td) ISvrmfcfc 


(3.11) 


only depends on three parameters: the two masses nikk, mx and the effective coupling 
Xa \yLN,a\'^■ Fixing the dark matter relic abundance to the observed best fit value for dark 
matter, = 0.1199, by Planck [43], we show in Fig. 1 two contour plots illustrating 

the dependence on the three parameters. We find that the effective coupling has to be of 
order 10“® for keV sterile neutrino masses in the range 2 — 100 keV and scalar doublets 
with electroweak-scale masses. Although there is no explanation for the smallness of the 
Yukawa coupling, it is technically natural. Larger scalar masses generally require smaller 
couplings to compensate for the suppression by the scalar mass rukk- Larger sterile neutrino 
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masses generally require larger effective couplings or larger scalar masses, because the DM 
abundance is proportional to the ratio mM/rukk- 


4 Free Streaming Horizon 


The free streaming horizon characterises the scale below which perturbations are suppressed 
in the power spectrum [44]. It is defined by the average distance a particle travels without 
any collisions 




rps = 


{v) 

a{t) 


dt 


(4.1) 


where {v) is the average velocity, tin denotes the time when the sterile neutrino is produced 
and to the time today. In order to evaluate it, we have to find the average velocity of the 
sterile neutrino 


w ^(1)^ 


1 for t < tpf^ 

1 


(4.2) 


where tnr {TN,nr) denotes the time (temperature of the sterile neutrinos) when neutrinos 
become non-relativistic, i.e. 

{p{TN,nr)) = rUN ■ (4.3) 

We can define the average momentum in terms of the momentum distribution function^ 
f{p,t) of the sterile neutrinos, which can be obtained from the Boltzmann equation 


m = c[f] 


(4.4) 


with the Liouville operator L[f] and the collision term C[f]. The Liouville operator is 
defined as [45] 

In analogy to the treatment in Refs. [15, 17], we introduce the dimensionless quantity 

PN 


XN = 


Tn ’ 


(4.6) 


which allows to write the distribution function as 

2\/9()7r ^ S'xhx 


fN{xN,TN) = 


E 


^ m\^JgZ{Td,X) Jo 


I 





dy\ —g(e 
XN 


,y+- 


with 


aiz) = - 


— --In 

l-hz 2 


z - 1 
z + 1 


(4.7) 


(4.8) 


The scalar particle mass (decay width) is denoted mx (hx) and we assumed that the 
effective number of relativistic degrees of freedom, g^, remains constant during production. 


^Note that we limit our discussion to a homogeneous isotropic Universe which is described by the 
Friedmann-Robertson-Walker metric and thus also take the distribution fnnction to be homogeneous. 
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Td,x is the decay temperature of particle X and we neglected the back reaction from 
inverse decays. See appendix App. B for a derivation of the distribution function. As the 
integrand is exponentially suppressed for y ^ 1, we can take the upper limit of the integral 
to infinity for temperatures Tjy <C mx/y/^xx- This is generally justified at late times, 
when pw r-u Tn < mx, thus we obtain for /^(xw) = limr^^o fN{xN,TN) 


fNi^x) 



X 


gxTxMpi 1 

'IXI’xV9* i'^d,x ) 






(4.9) 


In this approximation the function fx{x) has only one maximum at x ~ 1.54 and it falls 
off very quickly away from the maximum. A typical value for the momentum is thus of order 
xTtv. This justifies the approximation for small temperatures. The average momentum is 
given by 

(p(r^))~ 2.46 Tat (4.10) 

and thus lower compared to a sterile neutrino in thermal equilibrium with {p)t = S.ISTat. 
We find that the sterile neutrinos become non-relativistic at a temperature 


T'N,nr — 


mx 
2.46 ■ 


(4.11) 


The corresponding temperature of the SM thermal bath Tnr is obtained using the usual 
entropy dilution given in Eq. (C.15) with g^{Td) = g*{Td) = 110.75 and gl{Tnr) = 3.94 and 
thus the time tnr, when the sterile neutrinos become non-relativistic, is determined by 


t 


nr 


Ml ,, ‘-82 /siPk)j '''"Mu ^ 

V 9 *{Td) V 9i{Td) ) rn% \ mx J 


(4.12) 


Thus we can finally evaluate the integral for the free-streaming horizon by evaluating it 
piece-wise in the different regions set by tnr and teq, respectively 


^ ("s + In 0.047 Mpc (. (4.13) 

where feg = 1-9 x lO^^s and Ueq = 8.3 x 10“®. A detailed derivation can be found in 
Refs. [19, 20, 46]. Note that we do not have the additional entropy dilution factor, because 
it is already included in tnr in Eq. (4.12). Note, that the free-streaming horizon does 
not (strongly) depend on the mass of the heavy scalar or the effective coupling, which 
determine the sterile neutrino abundance. There is only an implicit dependence via the 
effective degrees of freedom at decay. A different scalar mass will lead to a different decay 
temperature, but for scalar masses above 100 GeV, the number of effective degrees of 
freedom stays almost constant. 

We show the free-streaming horizon in Fig. 2. The red region indicates when the 
free-streaming horizon becomes larger than 0.1 Mpc and the keV sterile neutrinos can 
be considered as hot dark matter. The blue region indicates the region when keV sterile 
neutrinos are in the cold dark matter regime. We take as benchmark value rps ^ 0.01 
Mpc following Refs. [19, 20]. 
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Figure 2: Free-streaming horizon vs. keV sterile neutrino mass. The regions of hot and 
cold dark matter are marked red and blue, respectively. A sterile neutrino mass niN = 7.1 
keV is indicated by the dashed orange line. 


We find free-streaming horizons within the desirable range for warm dark matter for 
a sizeable fraction of the parameter space. More precisely, the free-streaming horizon is 
in the warm dark matter range for keV sterile neutrino masses of 4 to about 53 keV. For 
example for mj\f = 7.1 keV, the keV sterile neutrino mass fitting the recently claimed X-ray 
line observation at 3.55 keV [37, 38], we obtain rps — 0.06 Mpc, which is well within the 
range of warm dark matter. 


5 Effective Degrees of Freedom 


With the introduction of an additional sterile neutrino, one might wonder about its con¬ 
tribution to the effective relativistic degrees of freedom. Its contribution can be quantified 


by[17] 


AiVeff = 


Ptv ~ TTij^n 
Plu 


with 


Piv — 


4 30 


(5.1) 


where we subtracted the non-relativistic energy density contained in the mass of the sterile 
neutrino. The energy density pM is given in Eq. (C.8). Thus we find that there are no 
additional relativistic degrees of freedom in the non-relativistic limit r]\f = m]\f/T]sf 3> x, 
but in the ultra-relativistic limit {r^ <C x) the additional effective relativistic degrees of 
freedom are 


AXeff(T) = 0.0379 


Mpi 

mkks/9* {Td) 



(5.2) 


Hence we find for a temperature T^, = Tbbn — 4 MeV shortly before the onset of big 
bang nucleosynthesis AAeff(rBBN) ~ 10“^ with S'KTbbn) = 10.75, gl{Td) = 110.75, 
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yLN,a ~ 10“^, the scalar mass rukk ~ 500 GeV, ~ 10 keV and hence a negligibly 
small contribution to the effective number of neutrinos. At recombination the temper¬ 
ature is well below the mass of the keV sterile neutrinos, misf, i.e. rjv 3> x, and thus 
^Aeff(rrec) — 0. 


6 Sterile Neutrino Mixing and the X-ray Line 

The scalar doublet scalar Hi, may obtain a vev, {Hi,) = ^0 , similar to the SM Higgs. 

A vev will induce mixing between the sterile neutrino and the active neutrinos 


yLN,aVv 

ruN 


( 6 . 1 ) 


and lead to a new contribution to the active neutrino mass via the seesaw mechanism 


^mv,ap = 1X1N sin 6a sin 9p . 


( 6 . 2 ) 


The mixing allows the production of the sterile neutrinos via neutrino oscillations [8-10], 
which can be estimated using the approximate formula [47] 


^N,osch^ ^ 0.2 X 


3 X 10-9 


rriN \ 1-8 

3kevJ 


(6.3) 


In addition the sterile neutrino can decay into a photon and an active neutrino generat¬ 
ing an X-ray line. This already constrains the mixing angle. See Ref. [11] for different 
constraints on a keV sterile neutrino which is produced via oscillations. Last year two 
independent groups observed an X-ray line at 3.55 keV [37, 38], which can be explained by 
the decay of sterile neutrino DM to a photon and a neutrino with an active-sterile mixing, 
sin^(20Q,) ~ 7 X 10“^^. However, the observation is still debated [48-52]. 

The vev Vi, can be naturally small and satisfy the X-ray bound, if it is induced via a 
small, possibly complex, Z 2 soft-breaking term [33, 36], 


Hsoft = + h.c.) 


(6.4) 


after the electroweak symmetry is broken. We obtain for the vev Vi, 

Re{nl2)v Im(Mi2) 


Re(uj, 


m: 


+ (-^3 + A 4 -|- A 5 ) 


Im('(;jy) ~ 


mn + (A 3 -I- A 4 - A 5 ) 


(6.5) 


There are no charge-breaking minima at leading order in /i^ 2 - The scalar masses will receive 
small corrections proportional to appearing at the second order in the scalar mass in 
Eq. (2.4). 

As the production via non-resonant neutrino oscillations can only give a subdominant 
contribution to the abundance of keV sterile neutrinos [11]^, the production is dominated by 
scalar decay. Thus we can translate the limit on the mixing angle from X-ray observations 

®This does not apply to sterile neutrino production via resonant neutrino oscillations [12]. Thus the 3.55 
keV X-ray line can be explained in models like in Ref. [53, 54]. 


- 9 - 









nikk [GeV] 

Figure 3: Contour plot showing the dark matter abundance Qh? for the benchmark point 
with niN = 7.1 keV and ^^sin^(20Q) ~ 7 x 10“^^. The red band indicates the 2(T-allowed 
region around the measured by Planck [43]. 


into a limit on the vev, Vjy. We find < 1 —10 MeV depending on the sterile neutrino mass 
rrijY ~ 2 — 100 keV and the scalar mass rrikk — 100 — 1000 GeV. Hence the contribution 
to the active neutrino masses in Eq. (6.2) is much smaller than the solar mass scale. The 
small splitting might explain the existence of pseudo-Dirac neutrinos, if the active neutrino 
mass originates from a Dirac mass term. Taking the claimed hint for an X-ray line at 
3.55 keV [37, 38], we show in Fig. 3 the dark matter abundance as a function of the scalar 
mass rukk and the effective Yukawa coupling ^/Yla \yLN,a\'^ fixing rriN = 7.1 keV and the 
active-sterile mixing ^^sin^(20a) ~ 7 x 10“^^. The red band indicates the 2iT-allowed 
region around the best-fit value measured by Planck [43]. The required effective Yukawa 
couplings are of order (5 — 22) x 10“® for scalar masses mkk — 100 — 1000 GeV. 

7 Conclusions 

In the Standard Model of particle physics extended by one singlet Majorana fermion, keV 
sterile neutrinos are usually produced via neutrino oscillations and the decay of the elec- 
troweak Higgs doublet scalar into the keV sterile neutrino is a subdominant contribution. 
However, we showed in this paper, that the decay of an electroweak doublet scalar can 
become the dominant production mechanism in a two Higgs doublet model, if the vev of 
the electroweak scalar doublet is small or even vanishes. Neutrino oscillations will only 
account for a minor additional contribution to the keV sterile neutrino abundance. Thus 
we do not consider a vev of the electroweak scalar doublet for the first part of the paper and 
only indicate in Sec. 6 the possible changes, when the electroweak doublet scalar obtains a 
vev. As long as the vev is small enough, the produced abundance can be approximated by 
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the result for a vanishing vev. The vev will lead to a mixing of the sterile neutrino with 
the active neutrinos. This renders it unstable via a two body decay allowing to search for 
the keV sterile neutrinos using X-ray line searches. 

We explicitly derived an analytic expression for the momentum distribution of the keV 
sterile neutrino for late times, studied its free-streaming horizon and briefly commented 
on its contribution to iVgg, the effective number of neutrinos in the early Universe, which 
is neglibly small. This production mechanism leads to a cooler spectrum of the sterile 
neutrino. The range for the sterile neutrino to be the warm dark matter is in between 4 
and 53 keV. 

The mechanism to produce sterile neutrinos via the decay of a Higgs doublet can 
be easily embedded in models of neutrino mass, which naturally explain the smallness of 
neutrino mass. This requires the addition of two or more sterile neutrinos. If the second 
Higgs doublet obtains a tiny vev, neutrino mass is naturally suppressed, while the Yukawa 
couplings of the other neutrinos can be relatively sizeable [33, 36]. For a vanishing vev of 
the second Higgs doublet, neutrino masses can be induced via the radiative seesaw in the 
scotogenic model [35] and dark matter is produced via freeze-in [34]. 

Note added: A recent publication [55] pointed out additional thermal corrections to 
the production rate, which become relevant at high temperatures. We do not expect 
any significant corrections, because keV sterile neutrinos are dominantly produced at low 
temperatures, T < mx- The inclusion of these effects goes beyond the scope of this work, 
which focused on a simple analytic discussion of this novel keV sterile neutrino production 
mechanism. A future numerical discussion should include these effects and also include the 
other neglected subdominant corrections discussed in Sec. 3. 
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A Decay Widths 


The decay width of the scalar fields are given as 


r(/c —)• NiUa) 

T{K^ ^ Niu^) 

r(A+ ^ NiQ 


mk\yLN,a? _ 
327r ^ 

'mKo\yLN,a\‘^ ^ 
327r ^ 

'mK+\yLN,a\‘^ ^ 
IGvr ^ 


_ mk\yLN,a? 
rn^ 327r 

_ ^^',2 _ rnxo\yLN,a\‘^ 
m^o 327r 


2 ^ rnK+\yLN,a\'^ 
IGvr 


(AT) 

(A.2) 

(A.3) 
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Hence the decay widths in the limit mw <C mx can be approximated by 

r(fc ^ NiVo,) _ r(K° ^ NiVcd _ T{K+ ^ NiQ) _ |j/Ljv,«P 

dk'f^k Sk'^^k^ SK+knx+ 32 tt 

B Dark Matter Momentum Distribution Function 


In case of the production of sterile neutrinos from decay of the real scalar k, the collision 
term is given by 

C[fN]^ ^ 


2Ex 


dUkdU,i2Tr)^6^^\pk - px - Pu)\M\ (A(l - /iv)(l - /.) - /7v/.(l + fk)) , 


with the Lorentz-invariant phase space element 

gd?p 


dll = 


(B.l) 


(B.2) 


(27r)32F ’ 

which is denoted dHy (dH^.) for neutrinos (the real scalar). The number of degrees of 
freedom are = 1 for the scalar and g^ = 2 for the neutrino. The distribution functions 
of the scalar (neutrino) are fk (fu) and the spin-averaged matrix element is given by 


_ 2 1 1 

\M\ = ^IvlM'^Pu ■ PN = {ml - mjf - ml) . 


(B.3) 


4'" ' " 4' 

Using the decay width T^ ~ \yLx\‘^mk/327r in the limit of negligible final state masses, 
mx = my = 0, we can rewrite the matrix element 


T-rvr2 SttT h 9 

\M\ = - -ml 

mk 

and generalise it to any of the three considered scalars using Eq. (A.4) 


(B.4) 


\M\ =STigxmxTx 


(B.5) 


Neglecting all terms proportional to fx and taking the ultra-relativistic approximation 
Ex — Pn, we obtain 

C[fx] = ‘^'^dxmxTx f ^nxdnj,(27r)"^(l(‘^)(px - Pn -Pu)fx (1 - fu) ■ (B.6) 

PN J 

Note that we do not neglect Pauli-blocking, i.e. do not approximate 1 — fy ~ 1. The 
integration over the neutrino momenta py is easily evaluated using the d-function of the 
momenta. As the integrand only depends on Ey and thus the scalar product px -PN, which 
can be evaluated via the delta function of the energies 


Ex - Ex = Ey = \py\ = \l'^x+d^N~ ■ PN , 

the collision term can we written as 


1 _ gxmxTx 
^[Jn\ — ■ ■ 


2p- 




N 


'PN + 


4pN 


T J 


Ex 

2T 


(B.7) 


(B.8) 
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Using the dimensionless variables 
PN,X 

XN,X = 


mx 


mx 


rj, , '^N = (B.9) 

we rewrite the Liouville operator in the radiation dominated epoch assuming to be 
constant 


Lifx] = Hr 


dfx 

dr 


(B.IO) 


and the collision term 


Cifx] = gx^x^^ j ^2 fxiux) (j-- fu dyx , (B.ll) 

where we use the Friedmann equation in a radiation dominated epoch to express the Hubble 
constant as H{T) = T'^/Mq = m\/Mor‘^ where Mq is defined in Eq. (3.4). Using our result 
from the collision term of a scalar decay, we hnd the Boltzmann equation for the distribution 
function fx of the sterile neutrinos 

dfx MogxTx 


dr 


2m\ 


J ,2 - u ■ (B.12) 


'^N •Jxn + 


Hence the distribution function can be obtained by a simple integration 

M„«rx r 


fx{xx,r) = 


2m\ 


.,r_l [ 

xlj Jx 


N 


\2 ) ) 


(B.13) 


If the scalar and the neutrino are in thermal equilibrium, they are described by a Bose- 
Einstein distribution function — 1 and Fermi-Dirac distribution fp}j{y) = 

+ 1, respectively. We can perform the integral over yx, 


fxixx.r) = 


Mpgx^x 

2m\ 




X 


N 


with 


9{z) = - 


1 


Changing the variables of integration to y = r'^/8xx we obtain 


- - In 
1 + z 2 

- M , 


z - 1 
z + 1 


fx{xx,r) = 


4V2Mogxrx 1 


m 


X 




y^gie 2 e^) dy 


We will employ the integral representation of the incomplete polylogarithm^ 


1 H 

/ 

r(s) Jb 


t 


s-1 


Lis{b, -z) = 


1 


i 


ed-jz - 1 

lS —1 


dt 




r(s) A eVz + l 


dt 


all z but Re(z) > 1 
all z but Re(z) < — 1 


(B.14) 

(B.15) 

(B.16) 

(B.18) 


^The integral representation of the incomplete poly logarithm is also known as the incomplete Bose- 
Einstein and Fermi-Dirac integrals, more precisely 


F,{x,b) = 


r(i 


J_ r 

i + 1 ) Jb 




dx = —Lij+i {b, —e^) . 


(B.17) 
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which are well-defined for Re(s) > 0, to rewrite the integral in terms of known functions. 
The usual polylogarithm is obtained for 5 = 0, i.e. Lis( 2 :) = Lis(0, z). The first summand of 
Eq. (B.15) can be integrated directly using the integral representation of the polylogarithm 
in Eq. (B.18), while the second term requires integration by parts to obtain fractions, which 
can then be integrated using Eq. (B.18). Finally the distribution function can be expressed 
at late times with r'^ ^ Sxjy by 


fNixN) = 


V^MogxTx 1 


m% 


2Li3 -e"^ This e"^ 


‘'X \/t/v - 

A comparison to the result without Pauli-blocking, 

V^MogxTx 1 


— Lis I —e 2 


M,nP / \ 

Jn i^N) = 


2m?^ 




Li, 


(B.19) 


(B.20) 


shows that the Pauli-blocking leads to the correction 


fNixx) - = 


2m\ 


1 


2Li3 I —e 2 I — Lis ( e 2 


— 3Lis ( —e 2 


(B.21) 


As expected the correction is negative and becomes larger for smaller xx- This agrees 
with the distribution function given in Ref. [56]. In the limit xx 3> 0 we obtain the 
Maxwell-Boltzmann result 




y/irMogx'^x e 


(B.22) 


which agrees with the result obtained from taking the Maxwell-Boltzmann limit in Eq. (B.13). 


fx^{xN,r) = 


^/xMogxTx e 


mx 

where we used the incomplete P function 


1 - 


3 \ 
2 ’ 4a; jv ) 

^W~ 


(B.23) 


poo 

r(s,x)= / . 

J X 


(B.24) 


The final expression for the distribution function at a finite temperature is given by 


fN{xx,r) = fx{xx) - 
with the function 


V^MogxTx 1 . / 


zh 


m^x \8xx 


,e 2 


(B.25) 


5,(5, z) = 2Li3 (5, z) + Lis (5, z) — Lis (5, —z) — _ 52 In 

2 2 2 Oyy TT 


4 ,3 , f - z\ 6=.oo 


e° + z 


0 . (B.26) 


We are mainly interested in the distribution function after the sterile neutrino DM produc¬ 
tion has finished and thus we will use the zero-temperature limit fx{xx) in the following. 
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The generalisation to multiple scalar fields is simply the sum of the different contribu¬ 
tions 


fNixN) 


^/90 


E 

X 


gxTxMpi 1 

9* {Td,x) 


2Li3 

2 



“h Lis 
2 




(B.27) 


The function has a maximum at x ~ 1.54 and falls off very quickly away from 

the maximum. A typical value for the momentum is thus of order x. After the sterile 
neutrinos are produced via the decay, they are completely decoupled from the thermal 
bath of the SM particles and generally have a different temperature Tx compared to the 
thermal bath of SM particles. Obviously the distribution function of the sterile neutrinos 
does not resemble a thermal Maxwell-Boltzmann, Fermi-Dirac or Bose-Einstein statistic, 
but represents a non-equilibrium distribution function. 


C Temperature of the Sterile Neutrino Dark Matter 


In the following we will calculate all relevant quantities for the sterile neutrino dark matter 
sector for sufficiently late time and low temperature T^r <C mx/^xx using the distribution 
function in Eq. (B.27). In order to simplify the notation we will also drop the subscript N 
from px and xx = Pn/Tx- The number density and the average momentum are given by 


n = gx J 

9\/5 


f ( \ '^N 


27r5/2 


(p) = ^ 

n 


- (|c( 4 )) E 


16 


x'^fx{x)dx 

gxTxMpi 


(C.l) 


X m\yJg'i{Td. 


=r3~ 0.091 

X) V 


gxTxMpi 


X 


Vg^d. 


--T] 


N 


X 


Pf.ip) = 


( 27 r 


2mr^ 


3. , ,29C(5)-32C(6)^ 

X jx{x)dx = —xxxnxxx 


'0 


26C(4)-32C(5)- 


2.46rjv , 
(C.2) 


where we used the integral formulas 

= a”'^^r(n -|- l)C(ra -|- s -|- 1) (C.3) 

x”LE(-e-^/“)dx = (-1 + 2-"-*) r(n + l)C(n + s + 1) (C.4) 

to evaluate the integrals. The result for the number density in Eq. (C.l) is consistent 
with the calculation of the yield Y = n/s in Eq. (3.7) using the integrated Boltzmann 
equation in the Maxwell-Boltzmann approximation. Using the expression for the number 
density and the distribution function we can determine how the temperature of the sterile 
neutrinos depends on the scale factor [57]. The sterile neutrinos are not interacting and 
thus the only effect is due to the cosmic expansion. After the abundance of the sterile 
neutrino is frozen-in, n{x)a? does not change and we can relate the number density n at 
time t to the number density n' at a later time t', n{x)a^ = n'{x'){a'Y. Then we obtain 
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from the number density of the sterile neutrinos with momenta between p and p + dp the 
following relation between the distribution functions at different times 


Tfjx^ffq{x)dx 



TSx'^fj^ix')dx' 


TnX^ f'N{^^x)dx 
a 1 : 


N 


(C.5) 


where we used the dependence of the momentum on the scale factor, p oc a to relate x 
and 

, _ P' _ a p Tn _ a Tn 
~ a'T mT'^~ a'T'/ 

Thus the distribution function f'jq at a later time t' has the same form with a temperature 

Tn = Tn^. (C.7) 

The expression for the energy density is similarly obtained from 


Pn = Pn 


^^^fN{p) = ^ y x‘^^Jx‘^ + r%fN{x)dx 


(p) n Tat <C X 
rriNn tat 3> x 


(C.8) 


with Tat = mi\f/Ti\f, where we took the ultra-relativistic or non-relativistic limit, respec¬ 
tively, in the last step. The pressure is given by 


Vn 



Tk 

Gvr^ 




fN{x)dx 


5 (p)n 

7 61C(6)-64C(7) PJV 

6 29C(5)-32C(6) 


Tn ^ X 
rw > X 


(C.9) 


and thus the equation of state is 



7 61C(6)-64C(7) 1 _ 1.16 
6 29C(5)-32C(6) 7^ “ 


Tat <C X 
tn ^ X 


(C.IO) 


showing that the sterile neutrinos behave like matter for temperatures much below rriN, 
i.e. Tat S> 1. Similarly the entropy density is given by 


SN 


PN + 'Pn 
Tn 


(1 -k w) 


PN 

Tn 


20 29C(5)-32C(6) 

3 26C(4)-32C(5) 

(1 + w) VnU 


n ~ 3.28n 


xat <C X 
Tn ^ X 


(C.ll) 


Thus we find that the entropy density scales like in the ultra-relativistic limit for 
xtv X and observe that the entropy of the sterile neutrino dark matter sector is separately 
conserved using Eq. (C.7). 

We can use the conservation of entropy to obtain a relation between the temperature of 
the SM thermal bath T and the temperature Tn of the sterile neutrino dark matter sector. 
Consequently we parameterize the entropy density similarly to the SM thermal bath 

27r^ q 27r^ q 

■s = -ntP^T and sn = -j^qI^nTn i (C.12) 


- 16 - 












we find for the entropy degrees of freedom in the sterile neutrino sector, gl 

s _ ■sp / 5 gxTxMpi 

9*,n ^ 50 Y gP(Td^j^) m\ 

Using the conservation of entropy 

gl{T,)Tlal = {gl{T)T^ + gl^^T%) 


(C.13) 


(C.14) 


and assuming that all scalars decay at the same time T^, we find for the temperature T/v, 


gl{T) \ ^ (g^\ 

T \9l{Td) - gl^N ) \gl{Td)) 


(C.15) 


In the non-relativistic limit, we can employ Eq. (C.7) to relate the temperature T^v to the 
temperature r/v,nr, 

TN=TN,nr^ (C.16) 

and ultimately the temperature of the SM thermal bath T. 
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